clear;
load('options_BE.mat')
load('Data_UK3_r.mat', 'Data_UK3_r')
nsect = length(Data_UK3_r);
%nsect=1
kjInd_UK3 = zeros(nsect,1);

for s = 1:nsect
    sigma=Data_UK3_r(s,3);    
    %sigma=1.0
    k_cool = @(k_j)sd_loge_int(k_j, sigma);
    x = fsolve(k_cool,5, options_BE)
    kjInd_UK3(s,1)=x;
    save('kjInd_UK3.mat', 'kjInd_UK3');
end



save('kjInd_UK3.mat', 'kjInd_UK3');


kappa1 = zeros(nsect,1);
kappa2 = zeros(nsect,1);
theta = zeros(nsect,1);

for s = 1:nsect
fun3 = @(z, k) (1-(z.^2)).*(z.*(exp(z))).^k .*exp(z);
Integral3 = integral(@(z)fun3(z, kjInd_UK3(s,1)),0,1);

kappa1(s,1)=kjInd_UK3(s,1).*exp(-kjInd_UK3(s,1)-1).*Integral3;
kappa2(s,1)=kappa1(s,1)./kjInd_UK3(s,1);
theta(s,1)=(kappa2(s,1).*(kjInd_UK3(s,1)+1))./((1./(kjInd_UK3(s,1)+1))-(kappa1(s,1)+kappa2(s,1))  );
end
thetajInd_UK3=theta;
save('thetajInd_UK3.mat', 'thetajInd_UK3');


init = ones(nsect,1);
CC=Data_UK3_r(:,4);
AA=theta;
ns=s;
beta_cool= @(B)eqsystem(AA, B, CC, ns);
    x = fsolve(beta_cool, init, options_BE)
    betajInd_UK3=x;
    save('betajInd_UK3.mat', 'betajInd_UK3');
    
    
intra_dist_UK3 = zeros(nsect,1);
intra_dist_UK3=100.*( ((kappa2.*((kjInd_UK3+1).^2)).^(-1./(kjInd_UK3+1))) -1  );
trick=betajInd_UK3'*thetajInd_UK3;
inter_dist_UK3 = zeros(nsect,1);
inter_dist_UK3=100.*((thetajInd_UK3./(trick))-1);

        save('inter_dist_UK3.mat', 'inter_dist_UK3');
        save('intra_dist_UK3.mat', 'intra_dist_UK3');